Analysis of ARPAV Temperature Data

library(tidyverse)
library(R2jags)
library(bayesplot)

global_seed <- 30172
set.seed(global_seed)

global_font <- "Roboto Condensed"
theme_set(theme_minimal(base_size = 15, base_family = global_font))
my_pal <- c("#EE9C39", "#20A47B", "#5A3920")

loc_path <- "./data/temp/malo"

temp_data <- NULL
for (type in c("min", "avg", "max")) {
    file <- list.files(
        path = loc_path,
        full.names = TRUE,
        pattern = paste0("*", type, ".csv")
    )

    tmp_df <- read.csv(file, sep = ";", na.strings = ">>")
    tmp_df <- tmp_df[-ncol(tmp_df)]
    names(tmp_df) <- c("year", 1:12)

    tmp_df <- tmp_df |>
        pivot_longer(-year, names_to = "month", values_to = "value") |>
        add_column(variable = type, .after = "month") |>
        mutate(
            date = make_date(year = year, month = month), .before = 1
        ) |>
        mutate(month = factor(
            month.abb[as.integer(month)], levels = month.abb)
        )

    temp_data <- rbind(temp_data, tmp_df)
}

temp_data <- arrange(temp_data, date) |>
    pivot_wider(names_from = "variable", values_from = "value") |>
    na.omit()

years <- unique(temp_data$year)

baseline <- temp_data |>
    filter(year <= years[5]) |>
    group_by(month) |>
    summarize(min = mean(min), avg = mean(avg), max = mean(max))

temp_data |>
    inner_join(baseline, by = "month") |>
    mutate(avg = avg.x - avg.y) |>
    ggplot() +
        geom_tile(aes(x = year, y = fct_rev(month), fill = avg)) +
        scale_x_continuous(
            breaks = scales::pretty_breaks(),
            expand = c(0, 0)
        ) +
        scale_fill_gradientn(
            colours = wesanderson::wes_palette("Zissou1", 10, "continuous")
        ) +
        labs(
            x = "Year", y = "Month", fill = "Deviation (°C)",
            title = paste0(
                "Deviation of the average monthly ",
                "temperature w.r.t. the ",
                years[1], "-", years[5], " period"
            )
        )

temp_data |>
    inner_join(baseline, by = "month") |>
    mutate(diff = (max.x - min.x) - (max.y - min.y)) |>
    ggplot() +
        geom_tile(aes(x = year, y = fct_rev(month), fill = diff)) +
        scale_x_continuous(
            breaks = scales::pretty_breaks(),
            expand = c(0, 0)
        ) +
        scale_fill_gradientn(
            colours = wesanderson::wes_palette("Zissou1", 20, "continuous")
        ) +
        labs(
            x = "Year", y = "Month", fill = "Deviation (°C)",
            title = paste0(
                "Deviation of the average monthly thermal ",
                "excursion w.r.t. the ",
                years[1], "-", years[5], " period"
            )
        )

smoothed <- temp_data |>
    group_by(year) |>
    summarize(temp = mean(avg))

# standard linear fit to get prior parameters
lm_fit <- lm(temp ~ year, data = smoothed)
lm_params <- summary(lm_fit)$coefficients[, 1]
lm_errs <- summary(lm_fit)$coefficients[, 2]

model <- sprintf("
    model {
        for (i in 1:length(year)) {
            mu[i] <- a + b * year[i]
            temp[i] ~ dnorm(mu[i], tau)
        }
        
        a ~ dnorm(%f, %f)
        b ~ dnorm(%f, %f)
        sigma ~ dexp(0.001)
        tau <- 1 / sigma^2
    }",
    lm_params[1], 1 / lm_errs[1]^2,
    lm_params[2], 1 / lm_errs[2]^2
)

pars <- c("a", "b", "sigma")

chain <- R2jags::jags(
    data = smoothed,
    inits = function() list(a = 0, b = 0, sigma = 1),
    parameters.to.save = pars,
    model.file = textConnection(model),
    n.chains = 3,
    n.iter = 100000,
    n.burnin = 2000,
    n.thin = 10,
    DIC = FALSE
) |> as.mcmc()

trace_plt <- mcmc_trace(
    chain, pars = pars,
    facet_args = list(nrow = 3, labeller = label_parsed)
) +
    facet_text(size = 15) +
    labs(title = "Traces") +
    bayesplot::theme_default(base_size = 15, base_family = global_font)

dens_plt <- mcmc_dens_overlay(
    chain, pars = pars,
    facet_args = list(nrow = 3, labeller = label_parsed)
) +
    facet_text(size = 15) +
    labs(title = "Densities") +
    bayesplot::theme_default(base_size = 15, base_family = global_font)

gridExtra::grid.arrange(trace_plt, dens_plt, ncol = 2)